Chapter 41
3D Surface Subdivision Methods

Table of Contents

41.1   Introduction

Subdivision methods are simple yet powerful ways to generate smooth surfaces from arbitrary polyhedral meshes. Unlike spline-based surfaces (e.g NURBS) or other numeric-based modeling techniques, users of subdivision methods do not need the mathematical knowledge of the subdivision methods. The natural intuition of the geometry suffices to control the subdivision methods.

Subdivision_method_3, designed to work on the class Polyhedron_3, aims to be easy to use and to extend. Subdivision_method_3 is not a class, but a namespace which contains four popular subdivision methods and their refinement functions. These include Catmull-Clark, Loop, Doo-Sabin and $$sqrt(3) subdivisions. Variations of these methods can be easily extended by substituting the geometry computation of the refinement host.

41.2   Subdivision Method

A subdivision method recursively refines a coarse mesh and generates an ever closer approximation to a smooth surface. The coarse mesh can have arbitrary shape, but it has to be a 2-manifold. In a 2-manifold, every interior point has a neighborhood homeomorphic to a 2D disk. Subdivision methods on non-manifolds have been developed, but are not considered in Subdivision_method_3. The chapter teaser shows the steps of Catmull-Clark subdivision on a CAD model. The coarse mesh is repeatedly refined by a quadrisection pattern, and new points are generated to approximate a smooth surface.

Many refinement patterns are used in practice. Subdivision_method_3 supports the four most popular patterns, and each of them is used by Catmull-Clark[CC78], Loop, Doo-Sabin and $$sqrt(3) subdivision (left to right in the figure). We name these patterns by their topological characteristics instead of the associated subdivision methods. PQQ indicates the Primal Quadtrateral Quadrisection. PTQ indicates the Primal Triangle Quadrisection. DQQ indicates the Dual Quadtrateral Quadrisection. $$sqrt(3) indicates the converging speed of the triangulation toward the subdivision surface.

The figure demonstrates these four refinement patterns on the 1-disk of a valence-5 vertex/facet. Refined meshes are shown below the source meshes. Points on the refined mesh are generated by averaging neighbor points on the source mesh. A graph, called stencil, determines the source neighborhood whose points contribute to the position of a refined point. A refinement pattern usually defines more than one stencil. For example, the PQQ refinement has a vertex-node stencil, which defines the 1-ring of an input vertex; an edge-node stencil, which defines the 1-ring of an input edge; and a facet-node stencil, which defines an input facet. The stencils of the PQQ refinement are shown in the following figure. The blue neighborhoods in the top row indicate the corresponding stencils of the refined nodes in red.

Stencils with weights are called geometry masks. A subdivision method defines a geometry mask for each stencil, and generates new points by averaging source points weighted by the mask. Geometry masks are carefully chosen to meet requirements of certain surface smoothness and shape quality. The geometry masks of Catmull-Clark subdivision are shown below.

The weights shown here are unnormalized, and $$n is the valence of the vertex. The generated point, in red, is computed by a summation of the weighted points. For example, a Catmull-Clark facet-node is computed by the summation of $$1/4 of each point on its stencil.

A stencil can have an unlimited number of geometry masks. For example, a facet-node of PQQ refinement may be computed by the summation of $$1/5 of each stencil node instead of $$1/4. Although it is legal in Subdivision_method_3 to have any kind of geometry mask, the result surfaces may be odd, not smooth, or not even exist. [WW02] explains the details on designing masks for a quality subdivision surface.

41.3   A Quick Example: Catmull-Clark Subdivision

File: examples/Subdivision_method_3/CatmullClark_subdivision.cpp

#include <CGAL/Cartesian.h>
#include <CGAL/Subdivision_method_3.h>
#include <iostream>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/IO/Polyhedron_iostream.h>

typedef CGAL::Cartesian<double>            Kernel;
typedef CGAL::Polyhedron_3<Kernel>         Polyhedron;

using namespace std;
using namespace CGAL;

int main(int argc, char** argv) {
  if (argc != 2) {
    cout << "Usage: CatmullClark_subdivision d < filename" << endl;
    cout << "       d: the depth of the subdivision (0 < d < 10)" << endl;
    cout << "       filename: the input mesh (.off)" << endl;
    return 0;
  }

  int d = argv[1][0] - '0';

  Polyhedron P;
  cin >> P; // read the .off

  Subdivision_method_3::CatmullClark_subdivision(P,d);

  cout << P; // write the .off

  return 0;
}

This example demonstrates the use of the Catmull-Clark subdivision method on a Polyhedron_3. The polyhedron is restricted in the Cartesian space, where most subdivision applications are designed to work. There is only one line deserving a detailed explanation:

Subdivision_method_3::CatmullClark_subdivision(P,d);

This example shows how to subdivide a simple Polyhedron_3 with Subdivision_method_3. An application-defined polyhedron might use a specialized kernel and/or a specialized internal container. There are two major restrictions on the application-defined polyhedron to work with Subdivision_method_3.

• Point_3 is type-defined by the kernel. Without Point_3 and the associated operations being defined, Subdivision_method_3 can not know how to compute and store the new vertex points.
• The primitives (such as vertices, halfedges and facets) in the internal container are sequentially ordered (e.g.  std::vector and std::list). This implies that the iterators traverse the primitives in the order of their creations/insertions.

Section 41.5 gives detailed explanations on those two restrictions.

41.4   Catmull-Clark Subdivision

When a subdivision method is developed, a refinement pattern is chosen, and then a set of geometry masks are developed to position the new points. There are three key components to implement a subdivision method:

• a mesh data structure that can represent arbitrary 2-manifolds,
• a process that refines the mesh data structure,
• and the geometry masks that compute the new points.

E. Catmull and J. Clark picked the PQQ refinement for their subdivision method, and developed a set of geometry masks to generate (or more precisely, to approximate) the B-spline surface from the control mesh. Subdivision_method_3 provides a function that glues all three components of the Catmull-Clark subdivision method.

template <class Polyhedron_3, template <typename> class Mask>
void PQQ(Polyhedron_3& p, Mask<Polyhedron_3> mask, int depth)

Polyhedron_3 is a generic mesh data structure for arbitrary 2-manifolds. PQQ(), which refines the control mesh p, is a refinement host that uses a policy class Mask<Polyhedron_3> as part of it geometry computation. During the refinement, PQQ() computes and assigns new points by cooperating with the mask. To implement Catmull-Clark subdivision, Mask, the geometry policy, has to realize the geometry masks of Catmull-Clark subdivision. depth specifies the iterations of the refinement on the control mesh.

To implement the geometry masks, we need to know how a refinement host communicates with its geometry masks. The PQQ refinement defines three stencils, and hence three geometry masks are required for Catmull-Clark subdivision. The following class defines the interfaces of the stencils for the PQQ refinement.

template <class Polyhedron_3>
class PQQ_stencil_3 {
  void facet_node(Facet_handle facet, Point_3& pt);
  void edge_node(Halfedge_handle edge, Point_3& pt);
  void vertex_node(Vertex_handle vertex, Point_3& pt);
};

Each class function in PQQ_stencil_3 computes a new point based on the neighborhood of the primitive handle, and assigns the new point to Point_3& pt.

We realize each class function with the geometry masks of Catmull-Clark subdivision.

template <class Polyhedron_3>
class CatmullClark_mask_3 {
  void facet_node(Facet_handle facet, Point_3& pt) {
    Halfedge_around_facet_circulator hcir = facet->facet_begin();
    int n = 0;
    Point_3 p(0,0,0);
    do {
      p = p + (hcir->vertex()->point() - ORIGIN);
      ++n;
    } while (++hcir != facet->facet_begin());
    pt = ORIGIN + (p - ORIGIN)/FT(n);
  }
  void edge_node(Halfedge_handle edge, Point_3& pt) {
    Point_3 p1 = edge->vertex()->point();
    Point_3 p2 = edge->opposite()->vertex()->point();
    Point_3 f1, f2;
    facet_node(edge->facet(), f1);
    facet_node(edge->opposite()->facet(), f2);
    pt = Point_3((p1[0]+p2[0]+f1[0]+f2[0])/4,
                 (p1[1]+p2[1]+f1[1]+f2[1])/4,
                 (p1[2]+p2[2]+f1[2]+f2[2])/4 );
  }
  void vertex_node(Vertex_handle vertex, Point_3& pt) {
    Halfedge_around_vertex_circulator vcir = vertex->vertex_begin();
    int n = circulator_size(vcir);    

    FT Q[] = {0.0, 0.0, 0.0}, R[] = {0.0, 0.0, 0.0};
    Point_3& S = vertex->point();
    
    Point_3 q;
    for (int i = 0; i < n; i++, ++vcir) {
      Point_3& p2 = vcir->opposite()->vertex()->point();
      R[0] += (S[0]+p2[0])/2;
      R[1] += (S[1]+p2[1])/2;
      R[2] += (S[2]+p2[2])/2;
      facet_node(vcir->facet(), q);
      Q[0] += q[0];      
      Q[1] += q[1];      
      Q[2] += q[2];
    }
    R[0] /= n;    R[1] /= n;    R[2] /= n;
    Q[0] /= n;    Q[1] /= n;    Q[2] /= n;
      
    pt = Point_3((Q[0] + 2*R[0] + S[0]*(n-3))/n,
                 (Q[1] + 2*R[1] + S[1]*(n-3))/n,
                 (Q[2] + 2*R[2] + S[2]*(n-3))/n );
  }
};

This example shows the default implementation of Catmull-Clark masks in Subdivision_method_3. This default implementation assumes the types (such as Point_3 and Facet_handle) are defined within Polyhedron_3. CatmullClark_mask_3 is designed to work on a Polyhedron_3 with the Cartesian kernel. You may need to rewrite the geometry computation to match the kernel geometry of your application.

To invoke the Catmull-Clark subdivision method, we call PQQ() with the Catmull-Clark masks we just defined.

PQQ(p, CatmullClark_mask_3<Polyhedron_3>(), depth);

Loop, Doo-Sabin and $$sqrt(3) subdivisions are implemented in the similar process: pick a refinement host and implement the geometry policy. The key of developing your own subdivision method is implementing the right combination of the refinement host and the geometry policy. It is explained in the next two sections.

41.5   Refinement Host

namespace Subdivision_method_3 {
  template <class Polyhedron_3, template <typename> class Mask>
  void PQQ(Polyhedron_3& p, Mask<Polyhedron_3> mask, int step);

  template <class Polyhedron_3, template <typename> class Mask>
  void PTQ(Polyhedron_3& p, Mask<Polyhedron_3> mask, int step);

  template <class Polyhedron_3, template <typename> class Mask>
  void DQQ(Polyhedron_3& p, Mask<Polyhedron_3> mask, int step)

  template <class Polyhedron_3, template <typename> class Mask>
  void Sqrt3(Polyhedron_3& p, Mask<Polyhedron_3> mask, int step)
}

The polyhedron class is a specialization of Polyhedron_3, and the mask is a policy class realizing the geometry masks of the subdivision method.

A refinement host refines the input polyhedron, maintains the stencils (i.e., the mapping between the control mesh and the refined mesh), and calls the geometry masks to compute the new points. In Subdivision_method_3, refinements are implemented as a sequence of connectivity operations (mainly Euler operations). The order of the connectivity operations plays a key role when maintaining stencils. By matching the order of the source submeshes to the refined vertices, no flag in the primitives is required to register the stencils. It avoids the data dependency of the refinement host on the polyhedron class. To make the ordering trick work, the polyhedron class must have a sequential container, such as a vector or a linked-list, as the internal storage. A sequential container guarantees that the iterators of the polyhedron always traverse the primitives in the order of their insertions. Non-sequential structures such as trees or maps do not provide the required ordering, and hence can not be used with Subdivision_method_3.

Although Subdivision_method_3 does not require flags to support the refinements and the stencils, it still needs to know how to compute and store the geometry data (i.e. the points). Subdivision_method_3 expects that the typename Point_3 is defined in the geometry kernel of the polyhedron (i.e. the Polyhedron_3::Traits::Kernel). A point of the type Point_3 is returned by the geometry policy and is then assigned to the new vertex. The geometry policy is explained in next section.

Refinement hosts PQQ and DQQ work on a general polyhedron, and PTQ and Sqrt3 work on a triangulated polyhedron. The result of PTQ and Sqrt3 on a non-triangulated polyhedron is undefined. Subdivision_method_3 does not verify the precondition of the mesh characteristics before the refinement.

For details of the refinement implementation, interested users should refer to [SP05].

41.6   Geometry Policy

Each geometry mask receives a primitive handle (e.g. Halfedge_handle) of the control mesh, and returns a Point_3 to the subdivided vertex. The function collects the vertex neighbors of the primitive handle (i.e. nodes on the stencil), and computes the new point based on the neighbors and the mask (i.e. the stencil weights).

This figure shows the geometry masks of Catmull-Clark subdivision. The weights shown here are unnormalized, and $$n is the valence of the vertex. The new points are computed by the summation of the weighted points on their stencils. Following codes show an implementation of the geometry mask of the facet-node. The complete listing of a Catmull-Clark geometry policy is in the Section 41.4.

template <class Polyhedron_3>
class CatmullClark_mask_3 {
  void facet_node(Facet_handle facet, Point_3& pt) {
    Halfedge_around_facet_circulator hcir = facet->facet_begin();
    int n = 0;
    Point_3 p(0,0,0);
    do {
      p = p + (hcir->vertex()->point() - ORIGIN);
      ++n;
    } while (++hcir != facet->facet_begin());
    pt = ORIGIN + (p - ORIGIN)/FT(n);
  }
}

In this example, the computation is based on the assumption that the Point_3 is the CGAL::Point_3. It is an assumption, but not a restriction. You are allowed to use any point class as long as it is defined as the Point_3 in your polyhedron. You may need to modify the geometry policy to support the computation and the assignment of the specialized point. This extension is not unusual in graphics applications. For example, you might want to subdivide the texture coordinates for your subdivision surface.

The refinement host of Catmull-Clark subdivision requires three geometry masks for polyhedrons without open boundaries: a vertex-node mask, an edge-node mask, and a facet-node mask. To support polyhedrons with boundaries, a border-node mask is also required. The border-node mask for Catmull-Clark subdivision is listed below, where ept returns the new point splitting edge and vpt returns the new point on the vertex pointed by edge.

  void border_node(Halfedge_handle edge, Point_3& ept, Point_3& vpt) {
    Point_3& ep1 = edge->vertex()->point();
    Point_3& ep2 = edge->opposite()->vertex()->point();
    ept = Point_3((ep1[0]+ep2[0])/2, (ep1[1]+ep2[1])/2, (ep1[2]+ep2[2])/2);

    Halfedge_around_vertex_circulator vcir = edge->vertex_begin();
    Point_3& vp1  = vcir->opposite()->vertex()->point();
    Point_3& vp0  = vcir->vertex()->point();
    Point_3& vp_1 = (--vcir)->opposite()->vertex()->point();
    vpt = Point_3((vp_1[0] + 6*vp0[0] + vp1[0])/8,
                  (vp_1[1] + 6*vp0[1] + vp1[1])/8,
                  (vp_1[2] + 6*vp0[2] + vp1[2])/8 );
  }

The mask interfaces of all four refinement hosts are listed below. DQQ_stencil_3 and Sqrt3_stencil_3 do not have the border-node stencil because the refinement hosts of DQQ and $$sqrt(3) refinements do not support global boundaries in the current release. This might be changed in the future releases.

template <class Polyhedron_3>
class PQQ_stencil_3 {
  void facet_node(Facet_handle, Point_3&);
  void edge_node(Halfedge_handle, Point_3&);
  void vertex_node(Vertex_handle, Point_3&);

  void border_node(Halfedge_handle, Point_3&, Point_3&);
};

template <class Polyhedron_3>
class PTQ_stencil_3 {
  void edge_node(Halfedge_handle, Point_3&);
  void vertex_node(Vertex_handle, Point_3&);

  void border_node(Halfedge_handle, Point_3&, Point_&);
};

template <class Polyhedron_3>
class DQQ_stencil_3 {
public:
  void corner_node(Halfedge_handle edge, Point_3& pt);
};

template <class Polyhedron_3>
class Sqrt3_stencil_3 {
public:
  void vertex_node(Vertex_handle vertex, Point_3& pt);
};

The source codes of CatmullClark_mask_3, Loop_mask_3, DooSabin_mask_3, and Sqrt3_mask_3 are the best sources of learning these stencil interfaces.

41.7   The Four Subdivision Methods

namespace Subdivision_method_3 {
  template <class Polyhedron_3>
  void CatmullClark_subdivision(Polyhedron_3& p, int step = 1) {
    PQQ(p, CatmullClark_mask_3<Polyhedron_3>(), step);
  }

  template <class Polyhedron_3>
  void Loop_subdivision(Polyhedron_3& p, int step = 1) {
    PTQ(p, Loop_mask_3<Polyhedron_3>() , step);
  }

  template <class Polyhedron_3>
  void DooSabin_subdivision(Polyhedron_3& p, int step = 1) {
    DQQ(p, DooSabin_mask_3<Polyhedron_3>(), step);
  }

  template <class Polyhedron_3>
  void Sqrt3_subdivision(Polyhedron_3& p, int step = 1) {
    Sqrt3(p, Sqrt3_mask_3<Polyhedron_3>(), step);
  }
}

The following example demonstrates the use of Doo-Sabin subdivision on a polyhedral mesh.

File: examples/Subdivision_method_3/DooSabin_subdivision.cpp

#include <CGAL/Cartesian.h>
#include <CGAL/Subdivision_method_3.h>

#include <iostream>

#include <CGAL/Polyhedron_3.h>
#include <CGAL/IO/Polyhedron_iostream.h>

typedef CGAL::Cartesian<double>            Kernel;
typedef CGAL::Polyhedron_3<Kernel>         Polyhedron;

using namespace std;
using namespace CGAL;

int main(int argc, char **argv) {
  if (argc != 2) {
    cout << "Usage: DooSabin_subdivision d < filename" << endl;
    cout << "       d: the depth of the subdivision (0 < d < 10)" << endl;
    cout << "       filename: the input mesh (.off)" << endl;
    return 0;
  }

  int d = argv[1][0] - '0';

  Polyhedron P;
  cin >> P; // read the .off

  Subdivision_method_3::DooSabin_subdivision(P,d);

  cout << P; // write the .off

  return 0;
}

41.8   Other Subdivision Methods

File: examples/Subdivision_method_3/Customized_subdivision.cpp

#include <CGAL/Cartesian.h>
#include <CGAL/Subdivision_method_3.h>

#include <cstdio>
#include <iostream>

#include <CGAL/Polyhedron_3.h>
#include <CGAL/IO/Polyhedron_iostream.h>

typedef CGAL::Cartesian<double>            Kernel;
typedef CGAL::Polyhedron_3<Kernel>         Polyhedron;

using namespace std;
using namespace CGAL;

// ======================================================================
template <class Poly>
class WLoop_mask_3 {
  typedef Poly                                         Polyhedron;

  typedef typename Polyhedron::Vertex_iterator         Vertex_iterator;
  typedef typename Polyhedron::Halfedge_iterator       Halfedge_iterator;
  typedef typename Polyhedron::Facet_iterator          Facet_iterator;

  typedef typename Polyhedron::Halfedge_around_facet_circulator
                                            Halfedge_around_facet_circulator;
  typedef typename Polyhedron::Halfedge_around_vertex_circulator
                                            Halfedge_around_vertex_circulator;

  typedef typename Polyhedron::Traits                  Traits;
  typedef typename Traits::Kernel                      Kernel;

  typedef typename Kernel::FT                          FT;
  typedef typename Kernel::Point_3                     Point;
  typedef typename Kernel::Vector_3                    Vector;

public:
  void edge_node(Halfedge_iterator eitr, Point& pt) {
    Point& p1 = eitr->vertex()->point();
    Point& p2 = eitr->opposite()->vertex()->point();
    Point& f1 = eitr->next()->vertex()->point();
    Point& f2 = eitr->opposite()->next()->vertex()->point();

    pt = Point((3*(p1[0]+p2[0])+f1[0]+f2[0])/8,
	       (3*(p1[1]+p2[1])+f1[1]+f2[1])/8,
	       (3*(p1[2]+p2[2])+f1[2]+f2[2])/8 );
  }
  void vertex_node(Vertex_iterator vitr, Point& pt) {
    double R[] = {0.0, 0.0, 0.0};
    Point& S = vitr->point();

    Halfedge_around_vertex_circulator vcir = vitr->vertex_begin();
    int n = circulator_size(vcir);
    for (int i = 0; i < n; i++, ++vcir) {
      Point& p = vcir->opposite()->vertex()->point();
      R[0] += p[0]; 	R[1] += p[1]; 	R[2] += p[2];
    }
    if (n == 6) {
      pt = Point((10*S[0]+R[0])/16, (10*S[1]+R[1])/16, (10*S[2]+R[2])/16);
    } else if (n == 3) {
      double B = (5.0/8.0 - std::sqrt(3+2*std::cos(6.283/n))/64.0)/n;
      double A = 1-n*B;
      pt = Point((A*S[0]+B*R[0]), (A*S[1]+B*R[1]), (A*S[2]+B*R[2]));
    } else {
      double B = 3.0/8.0/n;
      double A = 1-n*B;
      pt = Point((A*S[0]+B*R[0]), (A*S[1]+B*R[1]), (A*S[2]+B*R[2]));
    }
  }
  void border_node(Halfedge_iterator eitr, Point& ept, Point& vpt) {
    Point& ep1 = eitr->vertex()->point();
    Point& ep2 = eitr->opposite()->vertex()->point();
    ept = Point((ep1[0]+ep2[0])/2, (ep1[1]+ep2[1])/2, (ep1[2]+ep2[2])/2);

    Halfedge_around_vertex_circulator vcir = eitr->vertex_begin();
    Point& vp1  = vcir->opposite()->vertex()->point();
    Point& vp0  = vcir->vertex()->point();
    Point& vp_1 = (--vcir)->opposite()->vertex()->point();
    vpt = Point((vp_1[0] + 6*vp0[0] + vp1[0])/8,
                (vp_1[1] + 6*vp0[1] + vp1[1])/8,
                (vp_1[2] + 6*vp0[2] + vp1[2])/8 );
  }
};

int main(int argc, char **argv) {
  if (argc != 2) {
    cout << "Usage: Customized_subdivision d < filename" << endl;
    cout << "       d: the depth of the subdivision (0 < d < 10)" << endl;
    cout << "       filename: the input mesh (.off)" << endl;
    return 0;
  }

  int d = argv[1][0] - '0';

  Polyhedron P;
  cin >> P; // read the .off

  Subdivision_method_3::PTQ(P, WLoop_mask_3<Polyhedron>(), d);

  cout << P; // write the .off

  return 0;
}

The points generated by the geometry mask are semantically required to converge to a smooth surface. This is the requirement imposed by the theory of the subdivision surface. Subdivision_method_3 does not enforce this requirement, nor will it verify the smoothness of the subdivided mesh. Subdivision_method_3 guarantees the topological properties of the subdivided mesh. A genus-$$n 2-manifold is assured to be subdivided into a genus-$$n 2-manifold. But when specialized with ill-designed geometry masks, Subdivision_method_3 may generate a surface that is odd, not smooth, or not even exist.